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Abstract 

A simple dynamical cascade model for the evolution of free energy is considered in the context 
of gyrokinetic formalism. It is noted that the dynamics of free energy, that characterize plasma 
micro-turbulence in magnetic fusion devices, exhibit a clear predator prey character. Various key 
features of predatory prey dynamics such as the time delay between turbulence and large scale flow 
structures, or the intermittency of the dynamics are identifed in the quasi-steady state phase of 
the nonlinear gyrokinetic simulations. A novel prediction on the ratio of turbulence amplitudes 
in different parts of the wave-number domain that follows from this simple predator prey model 
is compared to a set of nonlinear simulation results and is observed to hold quite well in a large 
range of physical parameters. Detailed validation of the predator prey hypothesis using nonlinear 
gyrokinetics provides a very important input for the effort to apprehend plasma microturbulence, 
since the predator prey idea can be used as a very effective intuitive tool for understanding. 
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Predator-prey interactions are a very common paradigm in natural sciences, which pro- 
vide a powerful perspective for the interpretation of various complex phenomenajl|. In 
the context of fusion plasmas, the evolution of turbulence and self-regulating sheared flows 
that it drives, may eventually explain the dynamical coupling leading to the Low to High 
confinement (L-H) transition j^, yl in magnetic fusion devices. Interesting quasi-periodic 
activity, which may be linked to the predator-prey oscillations between turbulence and large 
scale flows in the form of zonal flows or geodesic acoustic modes, (GAMs) has been ob- 
served recently in a number of machines prior to and during the L-H transition J- 6|. The 
predator-prey dynamics also plays an important role in the nonlinear cascade process via the 
refraction of the turbulence in the low-fc (k being the wave-number) energy containing scales 
to high-Zc dissipative scales by the self-generated zonal flows j?-9|. This mediating role of 
zonal flows during the cascade of free energy has also recently been observed in gyrokinetic 
simulations 10j. 



Plasma turbulence in a strong magnetic field can be described by the gyrokinetic equation 
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-|l3|. which by filtering the rapid gyromotion, reduces the Vlasov equation to resolving 
a five- dimensional distribution function / = f(R,v\\,fj,,t), where R is the guiding center 
coordinate, v\\ is the velocity coordinate along the magnetic field Bq and \x = rriivj_/ '(2Bq) 
is the magnetic moment, which is an adiabatic invariant. In a field aligned geometry (R — > 
x,y,z with z the field aligned, x the radial and y the binormal coordinate), and separating 
between fluctuations and a Maxwellian equilibrium f = Fq + 5f, the gyrokinetic system of 
equations read: 



d t Sf k = iu* Ti F J $ k + iu D 5h k - D[Sh 
+ v Ti f ^ dzB ° d 5h k - v\\d z 5h k 



I 



]T (hxp.q) J %6h* (1) 

p+q=-k 

dpdv\\J 5h k = $ fc - , (2) 

where 6h k = 5f k + ^r 2 - Jo§ k has been used for compactness, the unknowns 5f k and are 
Fourier transformed in the plane (k x , k y ). Here, = (^fc)j where is the flux surface 
averaged electrostatic potential. Jo = Jo{k±p±) is the Bessel function of order zero, where 
k±(z) is the perpendicular wave number, and p±(z,fi) = Vti V ^f°(^) ~- Electrons are assumed 



adiabatic with temperature T e o and ions are characterized by an equilibrium density n^o and 



a temperature T^o (with T,o = T £ q in the following) , VTi = sj^ 2 -, — and Pi = fp- are 
respectively ion thermal velocity, cyclotron frequency and Larmor radius. Ion temperature 
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and density equilibrium profiles are contained in u*Ti = —jf- (l + \ V\\ + pB — |j j^- 
where: ^ = and ^- = -#k- . The drift frequency u D = _ V J^^1 ( Kk + Kk v ), 

contains magnetic unhomogeneity, details about the magnetic geometry can be found in 
Ref. D[5hk\ = c z d^5h k + c V]] d^^5h k + c±k 2 ™5h k corresponds to numerical dissipations. 

Simulations that we present in the following are performed using the GENE code Jj] 
Despite the fact that GENE is also adapted for electromagnetic and global problems 
the direct numerical simulations (DNS) presented here are restricted to local electrostatic 
ion temperature gradient driven turbulence (ITG) with adiabatic electrons. 

The gyrokinetic equation, as written in ([T|), has a number of nonlinearly conserved quan- 
tities, one of which is the so called free energy, which budget can be written as: 

d t £ = g-V, (3) 

where £ = n i0 j dA k ^^, Q = n i0 j dA k iw* Ti J Q $ k 8hl, and V = n i0 J dA k ^-D[5h k ] respec- 
tively define the free energy, its injection and dissipation (using the phase space integration 
f dA k = J2k x ,k v f ydzdv\\dp with the volume V = Y,k x ,k y f dz/B ). A local free energy 
balance in perpendicular Fourier space can be expressed as: 

d t £ i± =g e± +Af e± -V i± (4) 

where £± may be taken to correspond to any partition of the perpendicular Fourier space (for 
example k±£ < £± < k±e + i), and the contribution of the nonlinear term satisfies £Y J\ft ± = 0. 



In plasma turbulence, zonal flows |18| are of special importance, since these structures, 
extended over a given flux surface, play a regulating role on the turbulence that gener- 
ates them. The zonal flow free energy £ can be separated from the rest of the drift wave 
turbulence £ , where the energy budget takes the form: 

d t £ = A7-V, (5) 

d t £ = g +N -v. (6) 

It is important to note here that it is the free energy (which corresponds to "potential 
enstrophy" in the fluid limit) that is exchanged between the zonal flows and the drift waves 
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Figure 1. Time evolution of the normalized free energy £ associated with zonal flows and £ asso- 
ciated with the rest of the turbulence, for three different values of Ro/Lti = 6.0,6.92,8.0 from top 
to bottom. 



and not just the kinetic energy. The mechanism invoked here is not that of a classical 



"inverse cascade" but of a potential vorticity homogenization 



as disparate scale interactions in /c-space 
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which manifests itself 



2 if . Since there is no linear driving mechanism for 



the zonal flows (i.e. Q = 0), these structures feed on the free energy of the fluctuations and 
hence play the same regulating role on the underlying turbulence that a predator species 
plays on the population of a prey species. Figure [T] represents the time evolution of free 
energy (S and 8 have been normalized to their mean, and only a small fraction of the total 
time trace is represented in order to see the details of dynamics) , during the turbulent phase 
for three values of R n l Ln = 6.0, 6.92, 8.0, with other parameters being those of the ITG 



Cyclone Base Case 



22| (Ro/L n = 2.22, q = 1.4, s = 0.796, r = OT8.Ro, and T e0 = T,, 
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time delay in Rq/vtz 

Figure 2. Cross Correlation between Zonal Flow and drift wave turbulent signals, for different 
values of the logarithmic temperature gradient Rq/Lti- 

It can be observed that the free energy dynamics of the zonal flow and the turbulence are 
indeed largely correlated. In order to check if this dynamics exhibit predator-prey features, 
one can look at the phase relation between these two quantities to see if there exists a time 
shift between the turbulence and the zonal flow free energy fluctuations. In Figure [2} cross 
correlation between £ and £ is given as a function of the time lag. The average time delay 
At between £ and £ is given by the location of the maximal correlation in Figure [2j 

The predator prey type dynamics is also expected to have an intermittent nature. There- 
fore a look at the kurtosis is instructive: in Figure [3l the kurtosis associated with the free 
energy spectrum £ ky is plotted as a function of k y for the same runs as in Figure [TJ For 
all values of the temperature gradient, a clear separation is observed between the large and 
small scales, where the statistics of the small scales are very close to Gaussian (corresponding 
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to a kurtosis of 0), while the energy containing scales associated with wave vectors k y pi < 0.5 
show a clear departure from the Gaussian distribution, suggesting the presence of rare events 
with significant deviation from the mean. 

While the usual predator prey model already gives interesting perspective on the dynam- 
ics the fact that the predator-prey oscillations has the character of pulses in the fc-space 
transfer has to be considered also. Indeed, the free energy can be decomposed into the "en- 
ergy containing" component (corresponding to the energy injection scales), the "dissipative" 
component (corresponding to the dissipative scales) and a "Zonal Flow" component. These 
three domains are illustrated in Figure HI where the time averaged right hand side of the 
free energy equation i.e. Q ky — T> ky is given for various values of the logarithmic temperature 



gradient R^/L^i (the transfer term T ky is not represented since: T ky = T> ky — Q ky 23] ). 



hyperdiffusion c±k 4 Ll with c± = 0.375 
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Figure 4. Free energy spectrum £ ky for different values of the ion temperature gradient Rq/Lxi- 

Simulations have been performed on reduced grids (N x x N y = 48 x 24) by means of the 
GyroLES technique 24j] with a perpendicular hyperdiffusion model c± = 0.375. 

The spectral transfer character of the predator-prey dynamics ([5J [6]) can be studied by 
using a model of the form: 



d t E = AT — i/ F S 
<9 t £ 2 = A/" 2 - z^£ 2 



(7) 
(8) 
(9) 



where, £i is the free energy at the injection scale, £ 2 is the free energy at the dissipation scale 
and £ is the zonal free energy, and we have used the definitions: v F — T>/£, 7 = — 



and v = (T> 2 - G 2 )l£ 2 . 

This three domain model is very similar to the model studied in Ref. 0, except here 
we use free energy. As recently shown by Nakata et al. the free energy contribution 
of the nonlinear term can be expressed as a symmetrized triad transfer function: Ak = 
J2 P C^Sh^ShpShq, where C^' q is an operator converting the modified distribution func- 
tion Shk into the electrostatic potential This allows to write: 

W=Xhh 1 h 2l (10) 
Af 1 = \ 1 hh 1 h 2 , (11) 
M 2 = X 2 hh 1 h 2 . (12) 



where hi, h 2 and h can be defined for instance using the partition hs = J J s £ {k±) d 2 k± 
where the domain of integration in /c-space is chosen to correspond to the injection, dissipa- 
tion and zonal regions respectively. 

Following Ref. 0, equations ([71 El E]) can be averaged over the turbulent phase, so that 
time derivatives can be cancelled since the gyrokinetic simulation reaches a quasi-stationary 
state as shown in Figure [TJ Coefficients Ai, A2, A, 7, v and vp are constants in time, and by 
eliminating the product (hh\h 2 ) in the averaged equations, the two following relations can 
be obtained: 

An* (13) 



(£) A 7 
(£ 2 ) X 2 u F 



(14) 



(£) Xu 

The parameters Ai, A2 and A correspond to geometrical prefactors depending on the choice 
of the fc-space partition and the C , £' q that link $ and h, and not on the physical parameters. 
In contrast, the linear growth rate corresponding to 7, the small scale dissipation term u, 
and the zonal flow drag Vp are defined through Q\ j2 i G, ^1,2 an d T>, and thus are dependent 
on the free energy spectrum itself. Therefore, in principle the net dependence of these 
coefficients on physical parameters such as the temperature gradient can be quite nontrivial. 
This point is investigated in Figure 0, where the three parameters 7, vp and v are represented 
as functions of the imposed R /Lpi. 

In Figure [51 GyroLES simulations with various R^/Lpi are considered, with perpendicu- 
lar hyper diffusion amplitude c± = 0.5. The small scale dissipation v (dotted lines) is found 
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Figure 5. Effective Zonal Flows dissipation [yp , the solid line), growth rate (7 , the dashed line) 
and small scale dissipation (i/, the dotted line) as functions of Rq/Lti, with hyperdiffusion GyroLES 
(cjl = 0.5). 



approximately constant, except for small values of the temperature gradient. 7 and surpris- 
ingly vp present a nontrivial dependence with R /Lpi, very similar to the heat flux structure 



found in other studies 
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In order to verify Eqns. [13]and[T3J , the averaged free energy ratios (Si) / {£) and (£ 2 ) / (£) 
are represented as functions of the ratios vpjv and vp/'Y respectively in Figure [6l Six series 
of GyroLES simulations are considered, varying the diffusion or hyperdiffusion amplitudes, 
the temperature gradient Ro/Lpi, as well as the magnetic shear, for a total of 40 nonlinear 
GyroLES simulations. The advantage of using LES for these simulations, (apart from the 
gain in speed) is that it provides an easy handle on the small scale dissipation and allows 
us to modify v independently in order to explore the parameter space easily. 



The curves seems to agree with the theoretically predicted ratios in Eqns. [13] and [14] . 
Results are however found to depart from the theory when the turbulence level is decreased, 
especially in the low shear case and for very high perpendicular dissipation amplitudes. The 
deviation of the ratio 8\/8 from a straight line, is rather small, while it is more pronounced 
for 82/ S. This could be due to the reduction of the size of the dissipative range (which affect 
mainly 82) by the use of the GyroLES description. 

We have performed a detailed characterization of the dynamics and the associated spectral 
transfer, using the evolution of free energy in gyrokinetic turbulence with a partition of the 
fc-space corresponding to the scales of free energy injection, free energy dissipation, and 
large scale flow structures. Our results show that the free energy exchange between these 
components exhibits, any of the well-known characteristic features of the predator-prey 
dynamics. We have shown that the predicted relation between the average amplitudes of 
these different components as given in Eqns. [13] and [H] hold reasonably well in gyrokinetic 
simulations. 
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Figure 6. Free energy ratios (£i/£ top, bottom) as functions of the ratio between the effective 

Zonal Flow dissipation, the effective growth rate and the effective small scale dissipation (see Eqns[T3l 
andQ3J). Blue circles and red crosses stand respectively for diffusion and hyperdiffusion GyroLES 
model scanned along cj_, green plus and squares represent two Rq/Lti scans, with respectively 
cj_ = 0.5 and c± = 0.375 chosen for hyperdiffusion GyroLES. Black diamonds and stars correspond 
to a magnetic shear scan s with respectively c± = 0.5 and c± = 0.375. 
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